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I.  OVERVIEW  AND  8UHHARY 

Fundamental  to  this  work  has  been  the  development  of  a 
continuum  formulation  that  can  accurately  account  for  the  effects 
of  interlaminar  shear  and  interlaminar  normal  stress  variation 
through-the-thickness  of  a  laminate.  Furthermore,  emphasis  has 
been  particularly  on  tapered-twisted  airfoil  geometries  which  can 
be  analytically  represented  as  an  assemblage  of  thin  to  moderately 
thick  finite  elements.  To  achieve  solution  efficiencies,  various 
plate/shell  type  elements  have  been  developed  in  this  work  as 
opposed  to  the  more  computationally  intensive  solid  type  elements. 

On  the  basis  of  these  requirements,  alternative  continuum 

formulations  have  been  considered  and  are  herein  denoted  as  the 

(i)  Higher  Order  Displacement,  <ii)  Modif ied-Kirchhof f  and  (iii) 
Hybrid  Stress  formulations,  respectively.  "Shear  deformable" 
elements,  based  on  the  former  two  formulations,  have  been 
incorporated  in  a  computer  code  and  tested  on  the  basis  of 

correlations  with  known  analytical,  numerical,  and  experimental 
solutions.  Numerous  tests  have  been  performed  for  static,  dynamic 
and  buckling  behavior  of  laminated  structures.  Both  plate  and 
"doubly-curved"  shell  elements  have  been  formulated  and 
successfully  tested.  Note  that  use  of  a  shell  element  is  an 

especially  efficient  approach  in  modelling  airfoil  geometries. 

Significant  efforts  have  also  been  devoted  to  developing  a 
suitable  large  displacement  formulation.  Due  to  the  requirement 
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that  interlaminar  stresses  be  accurately  represented,  a  total 
LS9rangian  formulation  has  been  utilised  and  based  upon  the 
complete  Green's  strain  tensor.  A  geometric  and 
large-displacement  stiffness  formulation,  based  upon  a  form  of  the 
nonlinear  strain-nodal  displacement  relationship,  has  been 
developed  and  numerically  implemented  for  one  of  the  developed 
"shear  deformable"  plate  elements. 

Since  emphasis  in  this  work  focused  on  the  development  of 
incremental  response  solutions,  Including  damage  effects,  the 
computational  approach  needed  to  have  the  capability  to  (i) 
predict  and  differentiate  between  relevant  failure  modes,  (ii) 
modify  constitutive  equations  appropriately  and  (iii)  perform 
equilibrium  iterations  to  assure  stress  redistribution  based  upon 
the  extent  of  damage.  Use  of  "piecewise  smooth"  failure  criteria 
based  on  various  types  of  damage  has  provided  a  good  basis  for 
incrementally  tracking  damage.  This  approach  has  been 
incorporated  in  the  computer  code  using  various  stress  criteria. 
Alternative  strain-based  and  energy-based  approaches  have  been 
considered  including  the  representation  of  damage  via  the  use  of 
internal  state  variables. 

Note  that  integration  for  an  element  is  performed  on  a 
layer-by-layer  basis  which  allows  for  damage  effects  to  be 
characterized  at  the  layer  level.  Herein  a  layer  refers  to  either 
a  lamina  or  to  a  subset  of  adjacent  laminae  having  equal  ply 
orientations.  It  is  noteworthy  that  variation  in  strain  energy 
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can  be  calculated  as  damage  accumulates  and  that  such  information 
can  provide  the  basis  for  evaluating  sub-laminate  buckling  and 
even  delamination  growth. 

Experimental  data  has  been  utilised  to  substantiate  both 
stiffness  reduction  and  damage  predictions.  The  data  is  in  the 
form  of  failure  strengths  for  laminates  and  material  moduli 
variations. 

Computational  efficiencies  have  been  achieved  in  the 
numerical  procedures  primarily  by  using  an  equilibrium  formulation 
to  obtain  the  transverse  stresses.  This  approach  minimises  the 
need  for  3-D  solid  elements.  Furthermore,  a  "reduced  basis" 
approach  for  computing  nonlinear  dynamic  response  has  been 
developed  and  partially  implemented  in  the  computer  code. 

In  summary,  the  utility  of  the  developed  "shear  deformable" 
elements,  and  in  particular  the  use  of  an  equilibrium  formulation 
to  obtain  the  transverse  stresses,  has  been  demonstrated  for 
laminated  composite  geometries.  Realistic  structures  can  be 
modelled  with  the  use  of  a  mesh  generator,  which  has  been 
developed  in  order  to  generate  airfoil  geometries.  Since 
solutions  for  airfoil  response  were  not  available  with  which  to 
validate  the  computer  formulation,  laminated  composite  cylinders, 
with  various  boundary  conditions,  were  modelled  to  demonstrate  the 
capability  of  the  formulation  to  determine  transverse  stresses  in 
"curved"  geometries.  These  stresses  are  essential  in  tie 
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determination  of  damage  accumulation  because  of  the  typically  low 
interlaminar  strength  exhibited  by  laminated  composite  materials. 

Excellent  damage  predictions  have  been  obtained  by  coupling 
the  formulation  with  piecewise  continuous  stress  criteria  which 
provide  the  basis  for  differentiating  between  damage  modes.  These 
modes  include  fiber  fracture,  and  tensile-compressive  matrix 
cracking,  including  arrays  of  both  intralaminar  and  interlaminar 
cracks.  Strain  energy  calculations  have  been  performed  to 
determine  the  variation  in  strain  energy  release  rate  along  the 
boundary  of  an  interior  debonded  region.  These  calculations  are 
based  on  the  assumption  of  a  relatively  thin  sublaminatw  bounded 
by  a  rigid  "parent". 

The  numerical  formulation  developed  in  this  work  is 
especially  suited  to  analyzing  the  nonlinear  effects  and  damage 
progression  in  actual  structures.,  e.g.,  the  response  of  an 
airfoil  to  foreign  object  damage  or  the  response  of  a  helically 
wound  cylinder  to  a  tool  drop.  It  is  noted  that  "scaled  down" 
analyses  have  been  performed  for  primarily  two  reasons:  (1)  the 
lack  of  data  with  which  to  cor relate  the  results  and  (2)  the 
solutions  to  date  have  been  obtained  on  a  mini-computer,  while 
timely  solutions  for  the  examples  cited  would  have  to  be  obtained 
on  a  mainframe.  With  regard  to  the  second  point,  it  would  be 
worthwhile  to  either  "vectorize"  the  computer  formulation  or  to 
convert  the  code  to  a  parallel  processing  format. 
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The  sublaminate  analysis  could  be  sad*  more  raaliatic  by 
considering  boundary  conditions  that  reflect  the  actual 
displacement/slope  conditions  at  the  delamination  front.  This 
approach  is  achievable  through  more  coaputational  effort,  focused 
on  accounting  for  the  generalised  forces/displacement  that  occur 
along  the  boundary  of  the  debonded  region.  These  boundary 
conditions  could  be  supplied  by  the  coaplete  (global)  finite 
eleaent  aodal  to  a  local  aodel  of  the  delaainated  region. 

Finally,  there  are  strain-based  and  energy-based  damage 
criteria  that  have  the  potential  for  providing  aore  accurate 
daaage  predictions  than  those  obtained  herein  using  stress 
criteria.  Implementing  such  criteria  in  the  numerical  aodel, 
however,  would  require  significant  experimental  effort  to 
determine  certain  phenomonological  constants  needed  to  relate 
constitutive  behavior  to  damage  progression. 
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ZZ.  SUHMAJtY  SY  TASK 

This  ••ction  present*  technical  highlights  of  the  research 
efforts  to  date  for  each  of  the  three  tasks.  Details  of  the 
analytical  formulation  are  presented  in  the  Appendices. 

ZZ.l.  TASK  It  Nonlinear  Displacement  Formulation  for  Composite 

Media 

II. 1*1  Continuum  Formulation 

Two  variational  principles,  the  principle  of  minimum 
potential  energy  and  the  principle  of  modified  complementary 
energy,  are  generally  used  to  develop  two  distinctly  different 
finite  element  models,  the  assumed  displacement  model  and  the 
hybrid  stress  model  respectively.  These  models  incorporate  the 
effects  of  transverse  shear  and  normal  deformations  whose 
contributions  are  recognized  as  essential  for  accurate  laminate 
analysis  [1-10].  In  the  present  work,  emphasis  has  been  placed  on 
developing  displacement  based  models. 

Within  the  displacement  formulation,  element  stiffness 
matrices  are  determined  for  each  element,  these  matrices  are  then 
assembled  to  represent  the  final  system  of  equations  and  a 
solution  procedure  for  the  unknown  nodal  displacements  is 
provided.  Coordinate  transformations  to  describe  ply  orientations 
of  a  composite  media  are  taken  into  account.  The  in-plane 
stresses  are  calculated  from  constitutive  relations  of  orthotropic 
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continuum  whereas  transverse  shear  and  normal  atrataaa  art 
calculated  fro*  equilibrium  conside rations.  Finite  element  models 
have  been  tested  for  static,  dynamic  and  buckling  behavior.  The 
tast  problems  and  the  results  are  presented  in  Section  II. 1.4. 
The  finite  element  models  are  herein  briefly  discussed. 

A.  higher  Order  Displacement  Formulation 

The  through-the-thickness  effects  can  be  incorporated  into 
the  analysis  by  choosing  a  displacement  field  that  eliminates  two 
major  shortcomings  of  the  classical  plate  theory;  namely  normals 
remain  normal  and  in-plane  displacements  are  linear  through  the 
thickness.  These  shortcomings  are  eliminated  by  prescribing 
independently  the  reference  surface  displacements  and  rotations  of 
the  normal  and  including  higher  order  terms  for  in-plane 
displacements.  This  is  accomplished  in  the  plate  element 
formulation  by  the  following  variation 

2 

u(x,y, z )  -  uQ(x,y)  +  **x(x,y)  +  z  4x(x,y) 

2 

v(x,y, z )  -  vQ(x,y)  +  z*y(x,y)  +  z  fy(x,y) 
w(x,y,z)  -  wQ ( x , y ) 

The  neutral  surface  displacements  are  represented  by  uQ,vo  and  wo, 

the  rotation  about  the  y-axis  is  denoted  by  4>x  and  the  rotation 

2 

about  the  x-axis  is  ♦  .  The  coefficients  of  z  ,  i.e.,  *  and  ♦  , 

Y  x  y 

are  contributions  from  transverse  deformations  [5,6]  and  have  been 
found  to  be  significant  for  unsymmetric  laminate  geometries. 
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These  ttm  can  ba  omitted  in  tha  analytis  of  symmetric  laminates 
by  imposing  tha  appropriata  constraints. 

In  tha  shall  formulation,  tha  displacement  field  is  more 
complex  in  that  it  includes  components  of  tha  surface  unit-normal 
vector  and  three  rotational  components.  The  displacement 
components  become 

u(x,y,s)  -  uQ(x,y)  +  stNs*x<x,y)  -  Ny*#(x,y)J 
v(x,y,s)  -  vDU#y)  +  «(Nz*y(x,y)  +  Nx*#(x,y)l 
w(x,y,s)  -  wQ(x,y)  -  *lNx*x<x#Y>  ♦  Ny*y(x,y)) 

where  N  ,  N  and  N  are  components  of  the  surface  normal  at  a 

*  y  * 

particular  (x,y,z)  coordinate  location  for  the  element.  The  shell 

displacements  degenerate  to  those  of  a  plate  for  N  -  1  and  N  • 

2 

Ny  -  0  as  expected.  Note  that  Z  terms  are  not  included  in  the 

displacement  field  for  the  shell  element. 

The  elements  developed  are  designated  as  the  quadrilateral 
higher  order  displacement  (QHD)  models.  QHD40  is  an  eight-noded 
plate  element  with  seven  degrees  of  freedom  (three  midsurface 
displacements,  two  rotations  and  two  higher  order  terms  for 
in-plane  displacements)  per  corner  node  and  three  degrees  of 
freedom  (transverse  midsurface  displacement  and  two  rotations)  per 
mid-state  node.  QHD48  and  QHD48S  are  eight-noded  plate  and  shell 
elements  respectively,  with  six  degrees  of  freedom  (three 
midsurface  displacements  and  three  rotations)  per  node.  Element 
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QHD28  it  •  simplified  vtction  of  QHD40,  for  which  the  mid-side 
nodtt  art  eliminated.  It  should  be  noted  that  when  the  two  higher 
order  terms  for  in-plane  displacements  at  each  corner  node  are 
emitted,  QHD28  reduces  to  the  widely  used  four-noded  bilinear 
plate  eleaent  (QBD20) .  These  elements  produce  either  a  quadratic 
or  cubic  variation  in  the  transverse  shear  stress. 

B.  Aodif ied-Kirchhoff  Poraulation 

The  Kirchhoff-Love  assumption  for  normals  to  the  reference 
surface  is  relaxed  by  incorporating  shear  rotations  as  additional 
degrees  of  freedom  in  the  formulation  [10).  Thus  the  assumed 
displacement  field  allows  the  transverse  shear  deformations  but 
neglects  the  transverse  normal  deformations.  The  rotations  yx  and 
Yy  are  incorporated  in  the  displacement  variation  for  the  plate  as 
follows 


v(x,y)  -  wQ(x,y) 
u(x,y,z)  -  uQ(xfy)  -  z| —  +  yx 


raw 

v(x,y,z)  -  v  (x,y)  -  z —  «•  yv 

Uy  *) 


The  transverse  displacement  w(x,y)  is  chosen  such  that  it  will 
produce  stress  fields  that  characterize  the  transverse  effects 
accurately. 


This  approach  has  been  implemented  in  the  formulation  of  an 
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eight-node  quadrilateral  plate  element  with  32  degrees  of  freedom- 
QD32,  and  various  triangular  elements.  The  stress  fields  obtained 
for  these  elements  represent  a  quadratic  transverse  shear  stress 
variation.  The  quadrilateral  element  produced  the  best  results 
using  this  formulation. 


II. 1.2  Large  Displacement  Formulation 


Inclusion  of  geometrically  nonlinear  effects  in  the 
formulation  must  be  based  upon  the  geometry  to  be  analyzed  and 
upon  the  type  of  stress  prediction  capabilities  desired.  The 
classical  approach  to  thin  plate  analysis  has  been  to  use  the 
Kirchhof f-Love  assumptions  in  conjunction  with  the  nonlinear  von 
Karman  relations  [11,12].  As  previously  indicated,  the 
Kirchhof f-Love  assumptions  are  relaxed  in  this  work  to  allow  for  a 
more  accurate  definition  of  interlaminar  stress  variation.  These 
stresses  can  vary  substantially  through-the-thickness  for  the 
geometries  of  interest,  i.e.,  thin  to  moderately  thick  plate  type 
structures.  The  complete  Green's  strain  tensor  is  utilized  in 
this  work,  therefore,  to  account  for  all  significant  contributions 
to  the  interlaminar  stress  field.  With  respect  to  fixed  Cartesian 
coordinates,  x,  y,  and  z,  the  strain  tensor  has  the  form 


e 


x 


r 


xy 


dU 

-  + 

dX 

3u 

—  + 

ay 


■f3ul 

'3v' 

f*wl 

—  + 
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where  u,  v  and  v  represent  displacements  in  the  x,y,z  coordinate 
directions#  respectively.  Note  that  the  other  strain  components 
are  obtained  by  a  suitable  permutation.  In  small-displacement 
analysis#  the  quadratic  terms  are  neglected  to  give  simply  the 
linear  strain  approximation. 

Based  on  the  Green's  strain  tensor#  the  strain  to  nodal  point 
displacement  relationship  can  be  specified  for  elements  under 
development.  It  takes  the  form 

(t)  -  [BUM 

Where  { e >  is  the  vector  of  strain  components#  {A}  the  vector  of 
nodal  point  displacements  and  [B]  a  function  of  derivatives  of  the 
element  shape  functions.  The  quadratic  terms  in  the  strain  tensor 
result  in  [B]  being  a  function  of  displacement  state  and, 
therefore#  an  incremental  equilibrium  formulation  is  required. 
The  incremental  strain-nodal  displacement  relationship  takes  the 
form 


-  <  [BoJ  +  C  BL  J )  {5A} 

where  {5c}  and  {5A}  represent  incremental  strains  and  nodal 
displacements,  respectively,  [B0]  and  [BLJ  are  the  small  and  large 
displacement  contributions  to  the  incremental  strains.  Based  on 
the  incremental  equilibrium  equations,  the  displacement 
formulation  gives  the  force-displacement  relationships 
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tK-l  -  J  IB  r  fD] (B  ]  dv 

V 


|KL3  -  J  |tB0JT  [D][BL1  +  IBL]T[D][BL]  +Ibl]t[l][b0]  |dv 


where  [D]  is  an  elasticity  matrix  obtained  simply  from  the 
constitutive  equations  and  integration  is  over  the  volume  V  of  the 
element.  f K0 J  is  denoted  the  small-displacement  stiffness  matrix 
and  [kl]  is  denoted  the  large-displacement  stiffness  matrix. 
Since  response  is  also  a  function  of  stress  state,  the  geometrical 
stiffness  matrix  [K_]  is  required  and  is  obtained  from 


-  J  Ml 


IK  ]{5A)  -  J  MB  ]x{*)dV 
v 


where  {<r}  is  the  vector  of  stress  components. 


Inertial  effects  are  analytically  treated  as  a  mass  matrix 
[M]  which  is  a  function  of  density  and  the  element  shape 
functions.  These  matrix  forms  are  required  in  formulating 
static/dynamic  response  solutions  and  the  incremental  equilibrium 
equations  have  the  general  form 


[  M]  { $u}  +  ([KJ  +  IK-]  +  I  K_  ]  )  {Su}  -  { 6F) 


where  the  mass  and  stiffness  matrices  represent  an  assembly  of  the 
elemental  matrices  previously  discussed,  (Su)  and  (Su)  represent 
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the  incremental  displacements  and  accelerations  for  the 
mathematical  model  and  {8F}  represents  the  vector  of  incrementally 
applied  forces. 

A  complete  geometrically  nonlinear  formulation,  for  static 
and  dynamic  response  calculations,  has  been  implemented  in  the 
computer  code  only  for  the  QHD48  element.  Linear  elastic  buckling 
analysis  has  been  performed,  however,  with  each  of  the  elements, 
i.e.,  with  [*L]  omitted  from  the  equation  above. 

II. 1.3  Computer  Implementation 

Computer  coding  has  been  developed  for  the  purpose  of 
implementing  the  various  continuum  formulations.  The  code  has 
served  to  generate  static,  dynamic  and  buckling  solutions  using 
different  element  formulations.  All  of  the  element  integration  is 
performed  on  a'  layer-by-layer  basis  through-the-thickness  of  the 
laminate.  This  approach  is  fundamental  to  the  inclusion  of  damage 
mechanisms  in  the  formulation. 

Since  solution  of  the  equilibrium  equations  is  a  vital 

component  in  the  overall  solution  strategy,  it  is  appropriate  to 

discuss  the  numerical  methodology  used  in  solving  these  equations. 

The  intent  is  to  obtain  a  higher  ordered  variation  of  the 

transverse  shear  and  normal  stresses  (n  ,  a  ,  and  a  )  than  can 

xz  yz  zz 

be  obtained  via  the  constitutive  equations.  The  solution 
procedure  can  be  thought  of  as  described  below.  Assume  that  the 
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in-plane  stresses  (e  ,  e  ,  c  )  within  each  layer  of  a 

xx  yy  xy 

particular  element  have  been  determined  at  selected  locations, 
i.e.,  through  solution  of  the  constitutive  equations.  In  the  code 
as  presently  written,  these  locations  are  specified  as  the  element 
centroid  and  element  nodal  points.  The  equilibrium  equations  (in 
the  absence  of  body  forces)  have  the  indicial  form 

'ij.j  ■  0 

from  which  if  follows  that  the  through-the-thickness  shear  stress 
variation  can  be  written  in  numerical  form  for  the  ith  layer  as 

6°xzi  -  +  'xy.y’i  4Zi 

and 

Affyz .  ■  -( a  +  a  ).  AZ,  . 

J  i  xy,x  yy,y'i  i 

Here,  the  left-hand-side  represents  the  change  in  stress  from  the 
lower  to  the  upper  surface  of  the  ifc^'  layer  and  is  the 
thickness  of  the  ith  layer  at  a  particular  location.  The 
derivatives  with  respect  to  x  and  y  in  the  expressions  above  are 
readily  computed;  this  is  because  in-plane  stresses  within  a  layer 
are  related  to  element  displacements  through  derivatives  of 
element  shape  functions,  in  conjunction  with  a  material 
definition. 
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"*  -K  < 


,  L.f .  -_  A: 


For  an  n  layered  laminate,  n  aquations  can  ba  written  in 
tarms  of  both  the  unknown  shear  stresses  at  layer  interfaces  and 


the  shear  stresses  at  the  laminate  surfaces. 


Assuming  the 


laminate  has  shear-free  surfaces,  the  equations  above  give  n 
equations  in  n-1  unknowns,  so  that,  the  equation  set  is 
over-determined.  The  equations  have  the  matrix  form  below 


1 

-1  1 
-1  1 


n  x  (n-1)  (n-1)  x  1  (n  x  1) 


where  Ixzi  ■  “("xx  x  +  ffxy  y^iAZi  an<*  ffxzj  rePresents  the  shear 

i.  L  aL 

stress  acting  at  the  interface  of  the  j-1  and  j  layer.  A 
similar  equation  set  is  obtained  by  replacing  eX2j  with  «yzj  and 


Xxzi  with  1yzi * 


These  equations  are  solved  by  utilizing  a 


least-squares  orthonormalization  procedure  [13].  Due  to  the 
simplicity  of  the  terms  in  the  coefficient  matrix,  a  concise 
closed-form  solution  is  obtained.  Having  determined  the 
transverse  shear  stresses,  the  transverse  normal  stress  variation 
is  determined  through  the  numerical  form  of  the  third  equilibrium 
equation  of  the  ith  layer 


zz.  -  -(a  +a  ).6Z. 

i  '  xz,x  yz,y  i  i 


As  before,  the  left-hand-side  represents  the  change  in  stress 


through  the  ith  layer.  Appropriate  polynomial  functions  are 
utilised  to  describe  the  and  o„_  in-plane  variation.  These 

X  Z  jf  * 

functions  are  differentiated  to  obtain  the  right-hand-side  of  the 
equations  above.  Again  the  equation  set  is  overdetermined  because 
the  normal  tractions  are  known  at  the  laminate  surfaces.  Solving 
for  «  _  proceeds,  therefore,  in  identically  the  same  manner  as 
discussed  in  calculating  axz  and  ey2.  Parenthetically,  inclusion 
of  body  forces  can  be  accomplished  with  little  difficulty. 

As  a  final  note  with  regard  to  solving  the  equilibrium 
equations,  note  that  the  stresses  coming  from  the  constitutive 
equations  are  computed  in  an  "element”  coordinate  system.  In  the 
special  case  of  a  shell  element,  these  coordinates  are  generally 
not  directed  in  either  tangent-to  or  normal-to  the  shell 
directions.  An  additional  computational  requirement  in  this  case 
is,  therefore,  to  transform  the  "constitutive"  stresses  from 
element  to  shell  coordinates.  This  allows  numerical  integration 
to  proceed  on  a  layer-by-layer  basis  in  a  direction  normal  to  the 
shell  surface,  and  assures  that  the  stresses  are  determined  with 
accuracy. 

Successful  application  of  the  Higher  Order  Displacement  type 
elements  for  particularly  thin  geometries  requires  the  use  of 
reduced  numerical  integration.  This  approximation  technique 
brings  along  the  choice  of  implementing  it  overall  or  selectively 
to  the  strain  energy  components.  For  the  QHD  plate  formulations, 
only  the  transverse  shear  components  are  integrated  with  reduced 
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order  [14-16].  Good  results  are  obtained  with  the  shell 
formulations  by  under-integrating  all  strain  energy  terms. 

A  reduced  basis  numerical  algorithm  has  been  developed  and 
studied  with  respect  to  predicting  the  response  of  geometrically 
nonlinear  systems.  Good  results  have  been  obtained  for  some 
classical  problems.  The  computational  efficiency  of  this 
approach,  however,  is  only  realized  when  applied  to  a  very  large 
ordered  system  of  equations. 

A  preprocessor  has  been  developed  to  generate  the  finite 
element  mesh  for  an  airfoil  shape  of  multiple-circular-arc 
geometry.  Nodal  coordinates,  element  connectivity  and  unit 
surface  normals  are  all  generated  on  the  basis  of  limited 
geometrical  input  for  a  number  of  spanwise  locations  on  a  blade. 
The  normal  vector  definitions  are  needed  in  implementing  the 
QHD48S  shell  element.  Typical  mesh  geometries  are  shown  in 
Figures  1  and  2.  The  preprocessor  can  also  generate  finite 
element  meshes  for  the  special  case  of  a  cylindrical  shape,  which 
has  proven  very  useful  in  some  of  the  verification  testing 
presented  in  the  next  section. 

II. 1.4.  Analytical  Verification 

Significant  verification  results  relating  to  the  developed 
plate  and  shell  element  performance  have  been  presented  in 
previous  technical  reports  submitted  under  Contract  No. 
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F49620-82-K-0032.  Much  of  these  efforts  are  also  reported  in  the 
pipers  and  presentations  listed  in  Section  V  (Related  Activities) 
of  this  report.  Static,  dynamic  (fundamental  frequency  and 
transient)  and  buckling  response  predictions  for  various  laminate 
geometries  are  included  in  these  results  and  will  not  be  repeated 
herein.  It  is  noteworthy  that  the  QHD  formulation  gives  the  best 
results  overall  and  that  the  elements  are  suitable  for  the  study 
of  damage  accumulation  in  laminated  composite  structures. 

Results  comparing  the  performance  of  the  QHD4B  plate  and 
QHD48S  shell  elements  have  not  previously  been  presented.  Bach  of 
these  elements  has  six  degrees  of  freedom  per  node  and  is  similar, 
in  this  regard,  to  other  elements  found  in  the  literature.  The 
formulation  is  unique,  however,  due  to  the  layer-by-layer  approach 
to  obtaining  the  transverse  stresses  through  use  of  the 
equilibrium  equations.  This  is  the  common  thread  linking  all  of 
the  elements  developed  in  this  work.  For  the  QHD48  plate  element, 
3x3  Gaussian  quadrature  along  with  2x2  quadrature  for  the 
transverse  shear  components  is  employed.  The  element  exhibits  six 
rigid  body  modes  and  does  not  yield  any  undesirable  spurious 
modes.  In  the  shell  formulation,  2x2  Gaussian  quadrature  is  used 
to  integrate  all  strain  energy  terms.  This  numerical  approach 
produces  six  rigid  body  modes  plus  an  additional  zero-energy  mode. 
This  additional  mode  has  not  been  shown  to  have  a  deleterious 
effect  upon  element  performance. 

QHD48  and  QHD40  element  formulations  produce  essentially 
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identical  results  for  cylindrically  bent  and  rectangular  flat 
plates.  Xn  addition,  the  QHD48  element  has  been  used  to  aodel 
both  helically  wound  and  cross-plied  coaposite  cylinders.  First 
considered  is  a  long  cylinder  under  an  internal  pressure  of  50 
PSI,  with  inner  disaster  36  in.,  thickness  of  0.36  in.  and 
cylinder  length  of  »3  in.  The  laainate  geoaetry  is  given  as 
l35.3*,-35.3#,35.3#-35.3#,35.3e]g.  The  orthotropic  lamina 

properties  are  defined  as 

El  -  5. 136xl06  PSI 
Et  -  1 . 522xl06  PSI 
G,T  -  0.439xl06  PSI 
VLT  -  0.281 

and  typical  results  are  given  in  Figures  3  and  4.  It  is 
interesting  to  note  that  the  deformation  is  not  axisymmetric  as 
might  be  assuued.  Results  compare  quite  favorably  with  those 
obtained  using  a  3-D  solid  finite  element  formulation  [17]. 

The  QHD48  plate  and  QHD48S  shell  element  formulations  have 
been  utilized  to  calculate  the  transverse  shear  stress  variation 
that  occurs  in  the  'boundary  layer'  region  of  a  laminated 
cylindrical  shell.  Some  of  the  results  are  compared  to  those 
given  in  [18],  which  presented  an  analytical  solution  specifically 
for  determining  the  interlaminar  stresses  in  laminated  cylindrical 
shells.  Plate/shell  solutions  are  compared  to  those  given  in  [18] 
for  a  three-layer,  cross-ply  wound  cylinder.  Fiber  orientation  is 
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defined  at  |0*,  90*,  0*]  and  results  art  obtained  for  two 
different  sets  of  boundary  conditions  including  both 
slaply-supporttd  and  clasped  ends.  The  cylinder  geometry  and 
loading  is  the  saae  for  all  results  presented  herein.  The 
cylinder  has  a  length  L  -  50  inches,  radius  *  *  25  Inches  and  wall 
thickness  t  -0.25  inches.  Also,  the  cylinder  is  subjected  to  a 
unifora  internal  pressure  of  100  psi.  Two  aaterial  systeas  are 
considered  as  defined  below 

Boron-Epoxy: 


el 

■ 

32.5 

MPSI 

et 

m 

:  <84 

MPSI 

glt 

m 

0.642 

MPSI 

VLT 

- 

0.256 

Glaas-Epoxy: 

El  -6.00  MPSI 
Et  -  1.50  HPSX 
Glt  -  0.80  MPSI 
VLT  -  0.25  MPSI 

The  assumed  material  properties  and  geometry  definitions  are  all 
consistent  with  those  given  in  [18]. 

Since  the  laminated  geometry  under  consideration  is 
symmetric,  only  one  quadrant  of  the  cylinder  with  suitably 
specified  boundary  conditions  need  be  modelled.  Two  levels  of 
aesh  refinement  have  been  used  to  assure  that  convergent  solutions 
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h«v«  been  obtained.  These  include  both  10x10  end  16x16 
rectangular  meshes.  These  meshes  along  with  the  geometry  are 
shown  in  Figure  5.  Note  that  all  of  the  solutions  presented 
herein  are  for  the  aore  refined  aesh.  Calculated  values  for  the 
longitudinal  transverse  shear  stress,  ,  occurring  at  the 
interface  of  the  0*  and  90*  layers  are  given  in  Figures  6-9.  The 
first  two  of  these  Figures  relate  to  the  Boron-Epoxy  material 
systea,  while  the  latter  two  Figures  provide  results  for  the 
Glass-Epoxy  systea.  Regardless  of  which  boundary  condition  is 
considered,  the  results  indicate  essentially  aero  transverse 
stresses  away  froa  the  end  constraint.  Whereas  in  the  vicinity  of 
the  cylinder  end,  a  r boundary  layer'  develops  in  which  the 
transverse  stresses  becoae  quite  significant.  As  demonstrated  in 
these  Figures,  the  calculated  transverse  stresses  are  in  excellent 
agreeaent  with  the  analytical  solutions  given  by  Waltz  [18).  Note 
that  the  shell  element  model  gives  better  agreement  with  the 
analytical  solution  than  does  the  plate  element  model,  as  would  be 
expected.  While  the  10x10  aesh  results  are  not  presented  herein 
for  these  cases,  the  solutions  compare  well  with  the  16x16  results 
and  do  demonstrate  convergence. 

Finally,  calculated  transverse  stresses  are  presented  for  a 
cylindrical  shell  with  a  quasi-isotropic  layup  having  ply  sequence 
[0*,  +45",  -45°,  90° )g.  The  material  system  is  assumed  to  be  the 
Glass-Epoxy  previously  defined,  while  the  geometry  and  loading  are 
unchanged.  Figures  10  and  11  present  the  transverse  stress  'ryZ' 
variation  calculated  using  the  plate  and  shell  element 
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formulations,  again  for  both  clamped  and  simply-supported  boundary 
contfitlona.  Tha  'boundary  layer’  atraaaaa  are  similar  to  thosa 
shown  in  tha  pravious  cases.  Based  on  tha  previous  examples,  it 
is  believed  that  the  shell  element  provides  the  more  accurate 
results.  All  stress  results  represent  average  values  at  selected 
node  points  on  the  model.  it  is  interesting  to  note  that  the 
variation  in  these  stresses,  i.e.,  at  a  particular  node  shared  by 
more  than  one  element,  is  much  less  pronounced  in  the  shell  than 
in  the  plate  formulation.  It  seems  this  is  simply  a  further 
indication  that  the  shell  formulation  is  more  suitable  in 
modelling  structures  having  curvature. 

Some  verification  testing  of  a  developed  "reduced  basis" 
algorithm  has  also  been  performed  for  some  classical  problems. 
These  include  the  large  displacement  transient  response  of  both 
cantilevered  and  clamped-clamped  beams  subjected  to  step  loading. 
Results  are  given  in  Figures  12  and  13  demonstrating  the  accuracy 
obtained  in  these  problems  compared  to  directly  integrating  all  of 
the  dynamical  equations.  Note  that  the  basis  vectors  are 
comprised  of  both  Ritz  vectors  and  derivatives  of  Ritz  vectors. 
The  efficiencies  obtained  in  using  these  vectors  to  obtain 
solutions  for  linear  dynamical  systems  has  been  shown  in  [19,  20]. 
Reducing  the  order  of  the  equations  can  be  effective  when  the 
equations  do  not  need  to  be  updated  very  often  during  the  course 
of  a  transient  response  analysis.  Furthermore,  the  utility  in 
such  an  approach  is  only  realized  in  studying  the  response  of 
large-ordered  systems  of  equations. 
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12.2.  TASK  XX t  Incorporate  Dosage  Mechanisms  into  Dynaaic 

Keaponae  Formulation 

Relevant  failure  nodes  of  interest  include  those  listed  below 
(i)  fiber  fracture 
(ii)  fiber-matrix  debonding 

(ill)  matrix  cracking  (parallel  and  transverse  to  fibers) 
(iv)  delamination 

(v)  buckling  (possibly  at  layer  or  sub-laminate  level) 
Several  smooth  failure  criteria,  e.g.,  (21-24)  have  been  developed 
in  recent  years  to  represent  the  failure  of  composites.  These 
criteria,  to  varying  degrees,  can  predict  "failure"  but  do  not 
identify  a  particular  mode  of  failure.  In  performing  incremental 
"damage"  analysis,  it  is  essential  to  both  predict  failure  and  to 
characterize  it,  e.g.,  do  fibers  rupture,  does  delamination  occur, 
etc.  The  computational  approach  must,  therefore,  differentiate 
between  viable  failure  modes  and  appropriately  alter  the 

constitutive  equations  on  an  incremental  basis. 

Stress  Based  Approach 

One  approach  is  to  implement  a  piecewise  smooth  failure 
criteria,  e.g.,  (25)  in  the  finite  element  formulation.  The 
general  failure  criteria  is  then  comprised  of  m  separate 
inequalities  of  the  form 

Fj  ({*))<  1  J  j  -  1,2, ... ,m 


23 


at  the  layer  level  within  each  element.  These  criteria  should 
differentiate  between  (i)  tensile  and  compressive  fiber  failure, 
(ii)  tensile  and  compressive  matrix  failure  and  (iii)  delamination 
at  layer  interfaces  due  to  either  maximum  stress  or  buckling 
considerations . 

As  progressive  damage  occurs  throughout  incremental  loading 
(whether  it  be  static  or  dynamic),  it  is  essential  that  violation 
of  failure  criteria  inequalities  be  reflected  in  modification  of 
the  material  properties.  This  can  be  achieved  by  modifying  the 
appropriate  terms  in  the  constitutive  equations  to  reflect 
"stiffness  reduction".  When  the  strain  varies  between  tension  and 
compression,  as  in  the  case  of  transient  dynamic  response,  the 
numerical  model  must  reflect  the  differences  in  moduli  related  to 
whether  or  not  an  array  of  cracks  is  predicted  to  be  opened  or 
closed. 

In  conjunction  with  the  above  it  is  essential  to  perform 
equilibrium  iterations  within  each  analysis  increment.  This  is 
required  to  assure  that  stress  redistribution  is  properly 
accounted  for  as  damage  progresses. 

The  piecewise  smooth  failure  criteria  currently  implemented 
in  the  incremental  analysis  are  primarily  due  to  Hasbin  [25],  Lee 
[26],  Greszczuk  [27]  and  Hahn  [28].  Layer  stresses  are  defined  as 
shown  in  Figure  14  and  the  criteria  are  summarized  as  follows: 
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Fiber  Failure 


1.  Tension 

The  simplest  criterion  for  tensile  failure  of  a  composite  is 
the  maximum  stress  criteria.  The  failure  occurs  if  the  fiber 
tensile  stress  exceeds  the  allowable  normal  strength  eFN 

°L  *  *FN 


However,  this  is  a  drastic  approximation,  since  the  fibers  vary 
significantly  in  their  strength.  Lee  proposes  that  in  addition  to 
j  checking  this  criterion,  the  fibers  fail  if 

,2  .  2  .1/2  s 

°  LT  +  a  LZ  *  ffFS 

where  oFS  is  the  fiber  shear  strength.  On  the  other  hand,  the 
I  criterion  proposed  by  Hashin  for  the  tensile  fiber  failure  is 


°L  1 

“  1 

X  - 

/ 

2  + 

a2 

ffFN 

*  4 

T  1 

2 

a  LT  + 

\ 

°  LZ 

4 

2.  Compression 


For  compressive  loads  applied  along  the  fiber  direction,  a 
proposed  failure  mechanism  is  analogous  to  the  buckling  of  a 
column.  The  critical  fiber  buckling  stress  in  the  shear  mode  is 
given  by  Greszczuk  and  Hahn  as 

Gr 

a  *  - 

cs  ( 1-k ) 

where  Gf  is  the  resin  modulus,  and  k  is  the  volume  fraction  ratio. 
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Matrix  Failure 


1.  Tension 


The  composite  tensile  strength  transverse  to  the  fibers  is 
not  expected  to  deviate  significantly  from  the  matrix  tensile 
strength.  The  tensile  criteria  used  are  as  follows 


*•”  «T  i  °MN  0t  <”2TL  +  °2TZ)1/2  *  °( 


Hashin  j—  (a^+az)  +  j  “0T°Z^+  2  ^ ctLT  +  ffLZ 


where  (  oT  +  oz )  >  0 

The  matrix  normal  strength  and  the  matrix  shear  strength  are 
denoted  as  oM  and  crMC,  respectively. 


2.  Compression 


Under  compression,  failure  may  occur  by  shearing  along  a 
surface  through  the  matrix  parallel  to  the  fiber  axis.  The 
criterion  that  describes  this  is  proposed  by  Hashin. 


ffJ“l  - 


<aT+0  +- 


( °  TZ_ffT°Z) 


1 


1 


<  Vaz)+ 


1 


1 


where  9, 


is  the  compressive  matrix  strength. 


MNC 

Delamination 


Lee  proposes  that  delamination  occurs  if  either 

ffZ*ffDN  °r  <aLZ2  +  ffTZ2)1/2  -  ffDS 
where  «DN  and  <rDS  are  the  through-the-thickness  tensile  and  shear 

strengths  respectively.  Yet  another  form  used  by  the  authors  to 

identify  delamination  is  as  follows 


+ 


>  1 


Sublaninate  Buckling 


High  interlaminar  stresses  that  cause  delaminations  in 
composite  components  often  cause  localized  buckling  subsequent  to 
delamination.  These  high  stresses  may  promote  the  growth  of  the 
buckled  delaminated  region  and  lead  to  structural  failure  [29]. 
This  instability  related  crack  growth  can  be  studied  by  the 
virtual  crack  closure  technique  as  proposed  by  Whitcomb  [30]. 
Simply,  one  identifies  the  delaminaticn  zone  and  addresses  the 
layers  above  the  delamination  line  as  a  "sublaminate"  region  and 
the  layers  below  as  a  "parent”  region.  In  this  initial  model, 
parent  is  assumed  to  be  rigid.  Therefore,  one  models  only  the 
sublaminate  region  as  a  plate  with  clamped  boundaries.  Figure  15 
describes  the  geometry  of  the  model.  The  reactions  that  are 


calculated  at  the  boundaries  serve  as  forces  that  are  used  to 
close  the  debonded  region  (crack).  The  total  strain  energy 
release  rate  is  calculated  from  the  following  expression, 

G  -  0.5  •  (Nxex  ♦  Nyty  *  Nxytxy  +  MXKX  +MyFy  ♦  MxyKxy>. 

Then,  the  distribution  in  the  strain  energy  release  rate  is 
calculated  as  a  function  of  the  debond  geometry. 

Damage  Prediction  Calculations  (Stress  Based  Approach) 

The  damage  histories  for  selected  composite  laminates 
subjected  to  both  in-plane  and  bending  loads  have  been  determined. 
Note  that  both  static  and  transient  dynamic  loading  conditions 
have  been  considered. 

Uniaxial  Tension:  Response  to  Static  Load 

The  one  element  plate  model  of  Figure  16  is  employed  for  the 
uniaxial  tension  analysis  of  angle-ply  laminates.  The  laminate 
consists  of  three  layers  of  T300/5208  graphite  epoxy.  Two 
stacking  sequences  are  studied}  namely,  ie°/0o/-6°)s  and 
[ 0°/9O°/-e° ] s .  The  material  and  the  strength  properties  of 

T300/5208  are  given  in  Table  1.  The  first  and  last  ply  failure 
curves  of  the  two  laminates  as  a  function  of  Theta  are  shown  in 
Figures  17  and  18  for  a  statically  applied  load.  As  expected,  the 
first  and  the  last  ply  failures  for  uniaxial  laminates  of 
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(••/0 •/«•]«  of  Figure  17  occur  simultaneously,  except  for  angles 
greater  than  30*,  they  are  quite  separated.  The  ie*/90°/e*]s 
layup  of  figure  18  shows  that  for  9-60°,  75°,  and  90*  laminates, 
the  initial  and  final  failures  coincide,  whereas  for  angles  less 
than  60*,  they  are  easily  distinguished.  Table  2  displays  an 

alternate  view  of  the  damage  progression  in  a  (15°,  90*.  -15*1 
laminate,  where  initial  failure  occurred  at  the  60th  load 

increment  and  the  final  failure  at  the  117th  increment.  The 
column  headings  of  Table  2;  TF,  CF,  TM,  CM  and  DL  denote  Tensile 
Fiber  Failure,  Compressive  Fiber  Failure,  Tensile  Matrix  Failure, 
Compressive  Matrix  Failure,  and  Delamination,  respectively.  Note 
that  the  first  ply  failure  was  matrix  failure  in  the  90*  plys 

followed  by  the  fiber  failure  in  the  angle  plys.  Thus  one  can 
easily  identify  the  failure  mode  within  a  ply  for  a  given  load 

increment. 

Uniaxial  Tension:  Critically  Damped  Response  to  Step  Load 

The  same  uniaxial  tension  model  described  above  is  used  by 
the  incremental  dynamic  analysis.  To  verify  the  dynamic  response 
routines,  all  modes  in  the  transient  response  are  critically 

damped  to  approximate  an  incremental  static  solution.  The  first 
and  last  ply  failures  are  predicted  for  the  same  stacking 
sequences  as  described  in  the  incremental  static  analysis.  The 
results  for  the  [0®/O°/0®]  stacking  sequence  are  shown  in  Figure 
19.  The  first  and  last  ply  failures  occur  simultaneously  for 
0*0°,  15®  and  30*.  These  results  compare  quite  well  with  the 
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Tsai-Wu  failure  criteria  and  those  obtained  in  the  incremental 
static  analysis.  Likewise,  the  predicted  first  and  last  ply 
failures  for  the  stacking  sequence  [0°/90°/-9°]  compared 
exceptionally  well  with  the  Tsai-Wu  failure  criteria  and  with  data 
obtained  from  the  incremental  static  analysis  (Figure  20).  Table 
3  displays  the  progression  of  damage  within  the  l 15*/90#/-15# ] 

8 

laminate.  The  first  ply  failure  occurred  in  the  matrix  of  the  9C° 
plies  at  the  50th  load-time  increment.  The  final  failure  occurred 
in  the  fibers  of  the  15°  plies  at  the  97th  load-time  increment. 
The  damage  progression  predicted  by  the  incremental  dynamic 
analysis  is,  therefore,  quite  similar  to  that  obtained  using 
the  incremental  static  analysis. 

Uniaxial  Bar:  Transient  Response  to  Rectangular  Pulse 

A  [60°/0°/-60° lg  T300/5208  graphite/epoxy  bar  is  analyzed  for 
the  effects  of  damage  on  the  transient  response.  The  bar  has  a 
length  to  height  ratio  of  16.  Damping  is  not  included  in  this 
particular  analysis.  The  composite  bar  is  loaded  with  a 
rectangular  pulse  load  equivalent  to  45  percent  of  the  last  ply 
static  failure  load  for  the  laminate.  The  duration  of  the  pulse 
is  1.5  times  the  fundamental  period  of  the  bar.  The 
( 60°/0#/-60° ] s  composite  bar  is  analyzed  both  with  no  residual 
compressive  stiffness  and  with  ninety  percent  residual  stiffness 
after  a  tensile  failure  mode  has  occurred.  the  transient  response 
of  the  damaged  laminate  compared  to  the  linear  dynamic  response  is 
presented  in  Figure  21.  In  the  figure,  the  tip  displacement  is 
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normalised  by  the  maximum  amplitude  obtained  by  the  linear  dynamic 
aolution  and  time  is  normalised  by  the  fundamental  period  of  the 
bar.  At  point  A,  matrix  cracks  develop  in  the  60°  plies  in 
elements  1-7.  at  point  B,  the  maximum  tip  displacement  is  greater 
in  the  damaged  laminate  with  no  residual  compressive  stiffness 
than  in  the  damaged  laminate  with  ninety  percent  residual 
compressive  stiffness. 


Due  to  the  development  of  tensile  matrix  damage,  the 
amplitude  is  increased  in  the  first  cycle.  The  time  period  for 
♦  '»e  damaged  transient  response  also  tends  to  increase  as  damage 
lowers  the  stiffness  properties  of  the  laminate.  The  damaged 
laminate  with  no  residual  compressive  stiffness  exhibits  a  higher 
compressive  tip  displacement  than  either  the  linear  dynamic 
solution  or  the  damaged  laminate  with  ninety  percent  residual 
cot  essive  stiffness.  This  is  caused  by  the  lower  compressive 
sti  less  of  the  damaged  laminate. 


Tl  se  results  demonstrate  that  even  for  very  simple 
geometries,  it  is  essential  that  the  effective  material  moduli, 
due  to  the  opening/closing  of  crack  arrays,  be  realistically 
characterized  in  the  numerical  modelling  process. 


Four-Point  Bending:  Response  to  Static  Load 


The  bending  problem  of  Figure  22  is  modelled  with  four 
elements.  The  material  and  strength  properties  are  as  listed  in 
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Table  1.  The  laminate  is  unidirectional  and  consists  of 

twenty-four  layers.  For  the  bending  problem,  as  discussed  by 

Whitney  [31],  the  critical  aspect  ratio  is  defined  as  S-a  _  /x 

max  max  9 

and  for  the  present  geometry,  it  is  22.  Delamination  is  observed 
for  aspect  ratios  less  than  the  critical  in  the  middle  of  the 
laminate  as  the  load  is  increased.  Additional  matrix  and  fiber 
failure  accompany  delamination  as  shown  in  Table  4.  The 
interaction  curve  of  Figure  23  reveals  that  the  final  failure 
occurs  after  twenty  percent  load  increase  over  the  initial  failure 
load.  It  should  be  noted  that  for  aspect  ratios  less  than  the 
critical,  the  percent  increase  of  the  initial  failure  load  to  that 
of  final  failure  load  is  constant;  thus  if  one  reduces  the  shear 
strength  by  the  same  percentage,  the  final  failure  load  of  the 
corresponding  aspect  ratio  ends  up  on  the  interaction  curve.  This 
phenomenon  is  illustrated  by  the  dash-lines  of  Figure  11.  When 
the  aspect  ratios  are  higher  than  the  critical,  i.e.,  thin  plate, 
the  fiber  failure  at  the  outermost  laminae  proceeds  rapidly  toward 
the  center  and  within  four  percent  increase  of  the  initial  load, 
ultimate  laminate  failure  occurs.  A  typical  damage  progression  is 
displayed  in  Table  4  for  an  aspect  ratio  of  100.  It  is  also  of 
interest  to  observe  the  growth  of  failure  as  a  function  of  the 
total  strain  energy.  Figure  24  presents  the  variation  of  total 
strain  energy  in  each  lamina  of  a  [0°]^  laminate  with  an  aspect 
ratio  of  sixteen.  As  can  be  observed,  as  the  load  is  increased, 
the  strain  energy  increases  in  a  layer  until  that  layer  undergoes 
failure,  in  which  case  the  strain  energy  is  lower  than  its 
previous  value.  This  is  depicted  in  layers  12  -  20,  where  failure 


is  observed;  however,  the  undamaged  layers  1-11  and  21  -  24 
still  Maintain  higher  strain  energy  values  as  load  is  increased. 
Figure  25  focuses  on  the  saae  information  as  Figure  24,  yet 
provides  a  close-up  of  the  energy  variation  in  layers  5  through 
20.  As  discussed  above;  as  failure  is  detected,  the  strain  energy 
in  a  layer  decreases  and  causes  the  discontinuity  in  the  quadratic 
curve.  Note  that  the  failure  is  observed  as  starting  from  the 
Middle  layers.  On  the  other  hand,  for  an  aspect  ratio  of  100,  as 
shown  in  Figure  26,  the  failure  starts  from  the  outermost  laminae. 
In  this  case,  the  damaged  layers  are  1-4  and  20  -  24. 

Strain  Energy  Release  Rate  Calculations 

A  quarter  of  a  square,  debonded  sublaminate  region  is 
Modelled  a£  shown  in  Figure  27.  A  transverse  point  load  is 
applied  at  the  center.  The  strain  energy  release  rate  is 
calculated  from 

G  «  0.5*D11(M/P)**2 

where  is  the  bending  stiffness,  M  and  P  are  the  reaction 
moment  and  applied  transverse  load,  respectively. 

Simple  examples  are  chosen  to  evaluate  the  strain  energy 
release  rate  along  the  edge  of  the  debond,  for  isotropic  and 
orthotropic  layups,  namely;  [0°]4,  190°J4,  [45°  -45°jg.  Figure  28 
illustrates  the  distribution  of  G  along  the  y-axis  for  the  layups 
considered.  The  magnitude  of  the  total  strain  energy  is  the  area 
under  the  respective  curves  for  each  case.  For  all  cases,  it  is 
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obatrvtd  that  the  delamination  will  grow  at  tha  midsides  whara  G 
la  highast. 

■train  based  Approach 

An  altarnativa  to  modifying  tha  constitutiva  aquations  baaad 
on  violation  of  atrasa  failura  criteria,  is  to  use  a  atrain  baaad 
approach  that  reprasants  tha  accumulation  of  damage  through 
definition  of  internal  state  variables  in  conjunction  with 
experimentally  determined  phenomenological  constants.  The 
internal  state  variables  can  be  represented  as  vectors  having 
direction  (orientation  of  damage) and  magnitude  (extent  of  damage). 
This  approach  should  give  good  accuracy  in  characterising  the  type 
of  damage  that  develops  first,  i.e.,  intraply  matrix  cracking. 
This  approach,  combined  with  the  stress  criteria  for  fiber  failure 
and  delamination,  should  provide  more  accurate  predictions  than 
provided  by  use  of  the  stress  criteria  alone.  To  be  consistent 
with  the  finite  element  formulations  developed  in  this  work,  it  is 
important  that  damage  be  characterised  at  the  layer  level.  A 
laminate  damage  model,  in  combination  with  experimental  results 
for  certain  balanced  laminates,  can  serve  to  define  the  needed 
layer  damage  model.  These  models,  along  with  results  for  moduli 
variation  as  a  function  of  damage  accumulation,  have  been 
presented  and  discussed  in  the  previous  report  under  this  project. 
While  these  models  have  been  formulated,  they  have  not  been 
implemented  in  the  computer  formulation.  It  is  emphasized  that  a 
significant  amount  of  experimentation,  probably  using  non-standard 
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test  ipieium,  would  be  needed  to  fully  implement  such  an 
approach  to  modelling  daaagt  accumulation. 

11.3.3  TA8K  Ill:  Corralation  of  Poraulatad  Response  Model  with 

Experimental  Data 

Failure  results  presented  in  the  previous  section  do  compare 
well  with  the  experimental  data  given  in  (32).  Additional 
experimental  results,  while  not  plentiful,  do  exist  for  cases  in 
which  data  has  been  generated  to  quantify  the  effects  of  damage 
for  various  loading  conditions.  For  example,  the  extent  of  damage 
is  quantified  in  [33,34]  for  the  impact  response  of  composite 
pletes.  Correlation  of  the  present  formulation  with  such 

additional  data  would  be  a  worthwhile  endeavor. 
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Table  1,  Material  Proparti* a 


Elastic  Constanta 


Uniaxial  Tension 

Four-Point  Banding 

Ex  (GPa) 

138 

190 

E2  (GPa) 

10.6 

11 

E3  (GPa) 

10.6 

11 

G12  (GPa) 

6.4 

7.2 

G13  (GPa) 

6.4 

7.2 

®23  (GPa) 

6.4 

7.2 

V12 

0.3 

.38 

V13 

0.3 

.38 

V23 

0.3 

.38 

Strengths 

°FN  (MPa) 

1500  (1500)* 

1502  (1502) 

°FS  (MPa) 

68 

67.5 

(MPa) 

40  (246) 

41  (250) 

°hs  (MPa) 

68 

67.5 

(MPa) 

40 

41 

°DS  (MPa) 

68 

67.5 

*  Taras  in  parenthesis  ara  the  compressive  strength 


Table  2.  Damage  accumulation  of  laminate 

under  uniaxial  load  (static  analysis) 
115/90/-15), 


TF 

CF 

TM 

CM 

DL 

1 

117 

0 

0 

0 

0 

2 

0 

0 

60 

0 

0 

3 

117 

0 

0 

0 

0 

4 

117 

0 

0 

0 

.0 

5 

0 

0 

60 

0 

0 

6 

117 

0 

0 

0 

0 

3. 

Damage  accumulation 
under  uniaxial  load 
[ 15/90/-15)s 

of  laminate 
(dynamic  analysis) 

TF 

CF 

TM 

CM 

DL 

1 

97 

0 

0 

0 

0 

2 

0 

0 

50 

0 

0 

3 

97 

0 

0 

0 

0 

4 

o- 

0 

0 

0 

0 

5 

0 

0 

50 

0 

0 

6 

97 

0 

0 

0 

0 
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Figure  3 
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Figure  4 
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Fig.  5 .  The  reference  coordinate  system  and  the  two  meshes  including 
the  dimensions  of  the  dements  for  the  anisotropic, 
pressurized  cylinder. 


BORON  EP OXY,  CROSS-PLY,  PRESSURIZED  CYLINDER 
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Figure  6,  Transverse  shear  stress  at  the  interface  of  the  0  and  90  layers  for  a  boron-epoxy 
cylinder  with  simply-supported  edges* 


cylinder  with  clamped  edges 
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Figure  8.  Transverse  shear  stress  at  the  interface  of  the  0  and  90  layers  for  a  gla**s-epoxy 
cylinder  with  simply-supported  edges. 
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Figure  9.  Transverse  shear  stress  at  the  interface  of  the  0  and  90  layers  for 
cylinder  with  clamped  edges. 
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Figure  11.  Transverse  shear  stress  at  the  Interface  of  the  —45  and  90  layers  for  a  quasi-isotroplc 
glass-epoxy  cylinder  with  damped  edges. 
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Figure  w.  The  Stress  Components  in 
Natural  Coordinates 
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Figure  1 7.  Calculation  of  first  and  last  ply  failure 

for  a  [0/  0/  —  0], laminate-nonlinear  static  analysis 
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Figure  18.  Calculation  of  first  and  last  ply  failure 

for  a  [0/90/-  0], laminate-nonlinear  static  analysis 
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Figure  19 .  Calculation  of  first  and  last  ply  failure  for  a  critically 
damped  [5/0/  —  5]«laminate>nonlinear  dynamic  analysis 


Figure  20.  Calculation  of  first  and  last  ply  failure  for  a  critically 
damped  [5/90/  —  5Llaminate-nonlinear  dynamic  analysis 
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Figure  21  Effects  of  damage  on  the  transient  response  of  a  [60/0/  —  60] 
composite  bar  subjected  to  a  rectangular  pulse  load 


Figure  22.  Four-Point  Bending  Model 
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Figure  23.  Interaction  Curve  for  Four- Point  Bending  Load. 
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FIGURE  27.  Model  o£  Debonded  Sublamlnate  Region 


FIGURE  28.  G  Distribution  for  Sublaminate  Buckling 


